/**	@file Utilities.cpp
 *	@brief Auxiliary functions.
 */
#include "Utilities.h"
#include "Number.h"
#include "operators.h"
#include "Mul.h"
#include "Expr.h"
#include "ComparerIx.h"
#include "Substitution.h"
#include "Add.h"
#include "Pow.h"
#include "Fx.h"
#include "Sym.h"
#include <vector>
#include <algorithm>
namespace sympycpp {
Ex * Utilities::multiplierPair(Ex * e1, Ex * e2, const int level) {
	if ( e1->rsize() < e2->rsize()) {
		return multiplierPair(e2, e1, level);
	}
	if (e1->isAdd()) {
		if (e2->isAdd()) {// expand both
			Add * add1 = (Add *)e1;
			Add * add2 = (Add *)e2;
			Ex * tmp = new Number(0);
			add2->sign(add1->sign() ? add2->sign() : !add2->sign());
			for (iterator i = add1->args_.begin(); i != add1->args_.end(); ++i) {
				Ex * tmp1 = Operations::multiplication(*i, add2);
				tmp = Operations::addition(tmp1, tmp, STEALING);
			}
			expansion(tmp, level );
			return tmp;
		}
		else { //expand left
			Add * add = (Add *)e1;
			Ex * tmp = new Number(0);
			for (iterator i = add->args_.begin(); i != add->args_.end(); ++i) {
					Ex * tmp2 = Operations::multiplication(*i, e2);
					tmp = Operations::addition(tmp2, tmp, STEALING);
			}
			return tmp;
		}
	}
	else {
		if (e2->isAdd()) {//expand right
			Add * add = (Add *)e2;
			Ex * tmp = new Number(0);
			for (iterator i = add->args_.begin(); i != add->args_.end(); ++i) {
					Ex * tmp2 = Operations::multiplication(*i, e1);
					tmp = Operations::addition(tmp2, tmp, STEALING);
			}
			return tmp;

		}
		else {// do not expand
			return Operations::multiplication(e1, e2);
		}
	}
}
void Utilities::multiplier(Ex * & e, basic_container * vMul, const int level) {
	for (iterator i = vMul->begin(); i != vMul->end(); ++i) {
		Ex * tmp = multiplierPair(e, *i, level);
		delete e;
		delete *i;
		e = tmp;
	}
}
bool Utilities::expansion(Ex * & e, const int level) {
	if (level == 0) {
		return false;
	}
	if (e->isNum() || e->isSym()) {
		return false;
	}
	if (e->isFx()) {
		Fx * f = (Fx *)e;
		if (expansion(f->e_, level)) {
			e = f->create(f->e_,STEALING);
			e->sign(f->sign());
			f->e_ = 0;
			delete f;
			return true;
		}
		else {
			return false;
		}
	}
	if (e->isAdd()) {
		bool done = false;
		Add * addExp = (Add *)e;
		basic_container newAddends;
		for (iterator i = addExp->args_.begin(); i != addExp->args_.end(); ++i) {
			if (expansion(*i, level)) {// - 1)) 
				done = true;
				newAddends.push_back(*i);
				(*i) = 0;		
				if (!addExp->omit(i)) {
					Ex * tmp = 0;
					iterator j = addExp->args_.begin();
					if (i == j) {
						++j;
						expansion(*j, level);
						tmp = *j;
						(*j) = 0;
					}
					else {
						tmp = *j;
						(*j) = 0;
					}
					delete e;
					e = tmp;
					break;
				}
				--i;
			}
		}
		for (iterator i = newAddends.begin(); i != newAddends.end(); ++i) {
			e = Operations::addition(e, *i, STEALING);
		}
		return done;
	}
	if (e->isMul()) {
		bool done = false;
		Mul * mulExp = (Mul *)e;
		basic_container newMultiplicands;
		basic_container forMul;
		for (iterator i = mulExp->args_.begin(); i != mulExp->args_.end(); ++i) {
			if (expansion(*i, level - ((*i)->isPow() ? 0 : 1)) || (*i)->isAdd() ) {
				done = true;
				if ((*i)->isAdd()) {
					forMul.push_back(*i);
				}
				else {
					newMultiplicands.push_back(*i);
				}
				(*i) = 0;		
				if (!mulExp->omit(i)) {
					Ex * tmp = 0;
					iterator j = mulExp->args_.begin();
					if (i == j) {
						++j;
						expansion(*j, level - 1);
						tmp = *j;
						(*j) = 0;
					}
					else {
						tmp = *j;
						(*j) = 0;
					}
					tmp->sign(e->sign() == tmp->sign());
					delete e;
					e = tmp;
					break;
				}
				--i;
			}
		}
		for (iterator i = newMultiplicands.begin(); i != newMultiplicands.end(); ++i) {
			e = Operations::multiplication(e, *i, STEALING);
		}
		multiplier(e, &forMul, level);
		return done;
	}
	if (e->isPow()) {
		Pow * powE = (Pow *)e;
		const bool multTheorem = (powE->base_->isAdd() && powE->exponent_->isInteger());
		expansion(powE->base_, multTheorem ? level-1 : level);
		expansion(powE->exponent_, multTheorem ? level-1 : level);
		if (multTheorem) {
			Number * sum = (Number *)powE->exponent_;
			Sign s = plus;
			if (sum->lt(0)) {
				s = minus;
				sum->sign(plus);
			}
			Add * a = (Add *)powE->base_;
			Ex * newEx = multinomialTheorem(*a, *sum);
			delete e;
			if ( s == minus) {
				Ex * minOne = new Number(-1);
				e = Operations::power(newEx, minOne, STEALING);
			}
			else {
				e = newEx;
			}

		}
		return false;
	}
	return false;
}
void Utilities::makePermut(const int n, std::vector<int>  & sum, std::vector<std::vector<int> > & permut) {
	std::vector<int> v;
	for (int j = sum.size(); j < n; ++j) {
		sum.push_back(0);
	}
	sort(sum.begin(), sum.end());
	permut.push_back(std::vector<int>(sum));
	while(next_permutation(sum.begin(), sum.end())) { // STL
		permut.push_back(std::vector<int>(sum));
	}
}
void Utilities::nextSum(std::vector<int> & pre, std::vector<std::vector<int> > & res, const int n, const int sum) {
	if (n == 1) 
		return;
	for (int i = 1; i <= sum/2; ++i) {
		if(pre.size() == 0 || pre[pre.size()-1] <= i) {
			res.push_back(std::vector<int>(pre));
			res[res.size() - 1].push_back(sum - i);
			res[res.size() - 1].push_back(i);
			pre.push_back(i);
			nextSum(pre, res, n - 1, sum - i);
			pre.pop_back();
		}
	}
}
void Utilities::makeSums(const int n, const int sum, std::vector<std::vector<int> > & v) {
	std::vector<int> p;
	v.push_back(std::vector<int>());
	v[0].push_back(sum);
	nextSum(p, v, n, sum);
}
void Utilities::makeFactorials(const int x, std::vector<Number> & fs) {
	fs.push_back(1);
	for (int i = 1; i <= x; ++i) {
		fs.push_back(Number(1));
		fs[fs.size()-1].setToMultiplication(fs[fs.size()-2], i);
	}
}



void Utilities::multinomCoef(const std::vector<int> & sumands, const std::vector<Number> & f, Number & c) {
	for (size_t i = 0; i < sumands.size(); ++i) {
		c.setToDivision(c, f[sumands[i]]);
	}
}
void Utilities::adderHELP(const_iterator i, std::vector<Ex *> & src,  std::vector<std::pair<int, Ex*> > & help, const int index) {
	if ((*i)->isPow()) {
		Ex * exp = ((Pow *)*i)->exponent_->copy();
		help.push_back(std::make_pair(index, exp));
	}
	else {
		Ex * exp = new Number(1);
		help.push_back(std::make_pair(index, exp));
	}

}
void Utilities::adderSRC(const_iterator i, std::vector<Ex *> & src,  std::vector<std::pair<int, Ex*> > & help) {
	if ((*i)->isPow()) {
		Ex * base = ((Pow *)*i)->base_->copy();
		src.push_back(base);
		Ex * exp = ((Pow *)*i)->exponent_->copy();
		help.push_back(std::make_pair(src.size() - 1, exp));
	}
	else {
		Ex * expr = (*i)->copy();
		src.push_back(expr);
		Ex * exp = new Number(1);
		help.push_back(std::make_pair(src.size() - 1, exp));
	}
}
void Utilities::makeElements(const basic_container & p, std::vector<Ex *> & src, std::vector<std::vector<std::pair<int, Ex*> > > & helpS) {
	int iterer = 0;
	const_iterator pi = p.begin();
		
	if ((*pi)->isMul()) { // multiplicands contributes independent 
		const_iterator pp;
		if ((*pi)->isMultiple()) { // multiplicity is ignored
			pp = ((Mul *)*pi)->args_.begin();
			++pp;
		}
		else {
			pp = pi;
		}
		if ((*pp)->isMul()) {	
			for (const_iterator i = ((Mul *)*pp)->args_.begin(); i != (( Mul *)*pp)->args_.end(); ++i) {
				adderSRC(i, src, helpS[iterer]);
			}
		}
		else {
			adderSRC(pp, src, helpS[iterer]);
		}
	}
	else {
		adderSRC(pi, src, helpS[iterer]);
	}
	++iterer;
	++pi;
	while (pi != p.end()) { // transform all
		if ((*pi)->isMul()) { // multiplicands contributes independent
			const_iterator pp;			
			if ((*pi)->isMultiple()) {// multiplicity is ignored
				pp = ((Mul *)*pi)->args_.begin();
				++pp;
			}
			else {
				pp = pi;
			}
			if ((*pp)->isMul()) {
				for (const_iterator i = ((Mul *)*pp)->args_.begin(); i != ((Mul *)*pp)->args_.end(); ++i) {
					bool incompatible = true;
					int index = -1;
					for (std::vector<Ex *>::const_iterator j = src.begin(); (j != src.end()) && incompatible; ++j) {
						incompatible = INCOMPATIBLE == ComparerIx::compatibilityForMultiplication(*i, *j);
						++index;	
					}
					if (incompatible) {
						adderSRC(i, src, helpS[iterer]);
					}
					else {
						adderHELP(i, src, helpS[iterer], index);
					}
				}
			}
			else {
				bool incompatible = true;
				int index = -1;
				for (std::vector<Ex *>::const_iterator j = src.begin(); (j != src.end()) && incompatible; ++j) {
					incompatible = INCOMPATIBLE == ComparerIx::compatibilityForMultiplication(*pp, *j);	
					++index;
				}
				if (incompatible) {
					adderSRC(pp, src, helpS[iterer]);
				}
				else {
					adderHELP(pp, src, helpS[iterer], index);
				}
			}
		}
		else {
			bool incompatible = true;
			int index = -1;
			for (std::vector<Ex *>::const_iterator j = src.begin(); (j != src.end()) && incompatible; ++j) {
				incompatible = INCOMPATIBLE == ComparerIx::compatibilityForMultiplication(*pi, *j);	
				++index;
			}
			if (incompatible) {
				adderSRC(pi, src, helpS[iterer]);
			}
			else {// if both expressions are in multiplication, they will joined
				adderHELP(pi, src, helpS[iterer], index);
			}
		}
		++iterer;
		++pi;
	}
}
Ex * Utilities::makeAddend(const std::vector<int> & permutIndex, const std::vector<Ex *> & src, 
			   const std::vector<std::vector<std::pair<int, Ex*> > > & helpS, 
			   const Number & coef,
			   const std::vector<int> & indexMultiplicity,
			   const std::vector<int> & hasher) {
	Number * mult = new Number(1); 
	Number * pomMult = new Number(1);
	assert(permutIndex.size() == helpS.size());
	std::vector<Ex *> exponents(src.size(), 0);
	for (size_t i = 0; i < permutIndex.size(); ++i) {
		if (permutIndex[i] > 0) {
			for (size_t j = 0; j < helpS[i].size(); ++j) {
				int pos = helpS[i][j].first;
				if (exponents[pos] != 0) {
					Ex * expo =  new Number(permutIndex[i]);
					Ex * expo2 = helpS[i][j].second->copy();
					expo = Operations::multiplication(expo, expo2, STEALING);
					exponents[pos] = Operations::addition(exponents[pos], expo, STEALING);
				}
				else {
					Ex * expo =  new Number(permutIndex[i]);
					Ex * expo2 = helpS[i][j].second->copy();
					exponents[pos] = Operations::multiplication(expo, expo2, STEALING);
				}
			}
			pomMult->setValue(indexMultiplicity[i]);
			pomMult->setToPower(permutIndex[i]);
			mult->setToMultiplication(*mult, *pomMult);
		}
	}
	delete pomMult;
	basic_container multiplicands;
	if (exponents[hasher[0]] != 0) {
		Ex * base = src[0]->copy();
		Ex * addend = Operations::power(base, exponents[hasher[0]], STEALING);
		if (addend->isNum()) {
			((Number *)addend)->setToMultiplication(*(Number *)addend, coef);
			mult->setToMultiplication(*(Number *)addend, *mult);
			delete addend;
		}
		else {
			mult->setToMultiplication(*mult, coef);
			multiplicands.push_back(addend);
		}
	}
	else {
		mult->setToMultiplication(*mult, coef);
	}

	for (size_t i = 1; i < src.size(); ++i) {
		if (exponents[hasher[i]] != 0 ) {
			Ex * base = src[i]->copy();
			Ex * addend = Operations::power(base, exponents[hasher[i]], STEALING);
			multiplicands.push_back(addend);
		}
	}
	assert(multiplicands.size() >= 1);
	Ex * newExp;
	if (multiplicands.size() == 1) {
		newExp =  *multiplicands.begin();
	}
	if (multiplicands.size() > 1) {
		newExp = new Mul(multiplicands, STEALING);
	}

	if (! mult->eq(1)) {
		basic_container multiple;
		multiple.push_back(mult);
		multiple.push_back(newExp);
		return new Mul(multiple, STEALING);
	}
	else {
		delete mult;
		return newExp;
	}
	return 0;
}
Ex * Utilities::multinomialTheorem(const Add & base, const Number & exp) {
	const int N = base.args_.size();
	int intExp = exp.getInt();
	Number n(N);
	std::vector<std::vector<int> > sums;
	basic_container SUMA;
	std::vector<Number> factorials;
	
	makeSums(N, intExp, sums); //  
	std::vector<std::vector<std::pair<int, Ex *> > > helpS;
	std::vector<Ex *> src;
	for (size_t i = 0; i < base.args_.size(); ++i) {
		helpS.push_back(std::vector<std::pair<int, Ex *> >());
	}
	basic_container all_elements(base.args_);
	makeElements(all_elements, src, helpS);
	std::vector<Ex *> src2(src);
	std::vector<int> hasher(src.size());
	containerSort(src, ComparerIx::mullt);
	for (size_t i = 0; i < src.size(); ++i) {
		for (size_t j = 0; j < src.size(); ++j) {
			if (src[i]->str() == src2[j]->str()) {
				hasher[i] = j;
				src[i]->sign(plus);
				break;
			}
		}
	}
	std::vector<Sign> indexSign;
	std::vector<int> indexMultiplicity;
	for (const_iterator i = all_elements.begin(); i != all_elements.end(); ++i) {
		indexSign.push_back((*i)->sign());
		if ((*i)->isMultiple()) {
			indexMultiplicity.push_back(((Number *)*((Mul *)*i)->args_.begin())->getInt());
		}
		else {
			indexMultiplicity.push_back(1);
		}
	}
	makeFactorials(intExp, factorials);//
	for(size_t i = 0; i < sums.size(); ++i) {
		Number  coef(factorials[intExp]); // 
		multinomCoef(sums[i], factorials, coef); //DONE
		std::vector<std::vector<int> > permutIndex;
		makePermut(N, sums[i],permutIndex); // DONE
		for(size_t j = 0; j < permutIndex.size(); ++j) {
			Sign sign = plus;
			for (size_t i = 0; i < permutIndex[j].size(); ++i) {
				if (indexSign[i] == minus && (permutIndex[j][i] % 2) == 1) {
					sign = !sign;
				}
			}
			Ex * add1 = makeAddend(permutIndex[j], src, helpS, coef, indexMultiplicity, hasher);
			add1->sign(sign);
			SUMA.push_back(add1);
		}
	}
	containerSort(SUMA, ComparerIx::addlt);
	basic_container 					cleanSuma;
	iterator i = SUMA.begin();
	iterator j = SUMA.begin();
	++j;
	while(j != SUMA.end()) {
		if (ComparerIx::compatibilityForAddition(*i, *j) == INCOMPATIBLE) {
			cleanSuma.push_back(*i);
			i = j;
			++j;
		}
		else {
		Ex * helper = Operations::addition(*i, *j, STEALING);
			*i = helper;
			++j;

		}
	}
	cleanSuma.push_back(*i);
	Ex * result = new Add(cleanSuma, STEALING);
	return result;
}

void Utilities::taylor(const Ex * f, Ex * & fSeries, const Sym * byVar, const Ex * val, const int countDiff) {
	Ex * lastDiff = f->diff(*byVar); // f(x) 						
	Ex * byVar_val = Operations::subtraction((Ex *)byVar, val); //(x - a)			
	Number fact(1);
	Ex * substituted = f->copy(); //f(a)						
	Substitution::sub(substituted, (Ex *)byVar, val); 
	basic_container seriesX;
	seriesX.push_back(substituted);
	for (int i = 1; i <= countDiff; ++i) {
		Ex * nextDiff = lastDiff->diff(*byVar);				//	+
		substituted = lastDiff;						
		lastDiff = nextDiff;				
		Substitution::sub(substituted, (Ex *)byVar, val); // f'(a)
		Number exp(i);
		Ex * powBrackets = Operations::power(byVar_val, &exp); // (x-a)^i 	
		fact.setToMultiplication(fact,i);	// i!
		Ex * factorial = fact.copy();
		Ex * divSubFac = Operations::division(substituted, factorial, STEALING);
		Ex * member = Operations::multiplication(divSubFac, powBrackets, STEALING);
		seriesX.push_back(member);
	}
	delete lastDiff;
	delete byVar_val;
	fSeries =  constructorAdd(seriesX, STEALING);
}
Expr Utilities::set0_2Pi(const Expr & e1, bool & done) {
	Ex * newEx = e1.e_->copy();
	set0_2Pi(newEx, done);
	return Expr(newEx, Expr::STEALING);
}
void Utilities::set0_2Pi(Ex * & e_, bool & done) {
	done = false;
	if (e_->isMultiple()) { 
		iterator it = ((Mul *)e_)->args_.begin();		
		++it;
		if ((*it)->str() == pi.str()) {
			done = true;
			Number * multiplicity = (Number *) *((Mul *)e_)->args_.begin();
			multiplicity->sign(e_->sign());
			e_->sign(plus);
			const Number two(2);
			while ( two.lt(*multiplicity) || two.eq(*multiplicity)) {//p >= 2) 
				multiplicity->setToSubtraction(*multiplicity, two);
			}
			while  (multiplicity->lt(0)) {
				multiplicity->setToAddition(*multiplicity, two);
			}
			if (multiplicity->eq(0)) {
				(*((Mul *)e_)->args_.begin()) = 0;
				delete e_;
				e_ = multiplicity;
			}
		}
	}
	if (e_->isSym() && ((Sym *)e_)->title() == pi.str()) {
		done = true;
		e_->sign(plus);
	}
	if (e_->isNum() && ((Number *)e_)->eq(0)) {
		done = true;
	}
}

Ex * Utilities::constructorAdd(basic_container & ar, const allocationPolicy flag) {
	assert(ar.size() > 0);

	containerSort(ar, ComparerIx::addlt);
	
	basic_container cleanSuma;
	iterator i = ar.begin();
	iterator j = ar.begin();
	++j;
	bool nextZero = (*i)->isNum() && ((Number *)(*i))->eq(0);
	while(j != ar.end() && nextZero) {
		delete (*i);
		(*i) = 0;
		i = j;
		++j;
		nextZero = (*i)->isNum() && ((Number *)(*i))->eq(0);
	}	
	while(j != ar.end()) {
		if (ComparerIx::compatibilityForAddition(*i, *j) == INCOMPATIBLE) {
			cleanSuma.push_back(*i);
			*i = 0;
			i = j;
			++j;
		}
		else {
			*i = Operations::addition(*i, *j, STEALING);
			++j;

		}
	}
	cleanSuma.push_back(*i);
	*i = 0;
	i = cleanSuma.begin();
	if (cleanSuma.size() == 0) {
		return 0;
	}
	if (cleanSuma.size() == 1) {
		return *i;
	}
	return new Add(cleanSuma, STEALING);
}
}
